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INTRODUCTION 


Coordinate system selection is an important consideration in 
the time asymptotic numerical solution of any fluid flow or heat 
transfer problem. In solving such transient problems, the physical 
domain is usually transformed into a rectangular region with bound- 
aries coincident with the physical boundaries. Once this trans- 
formation is completed, the transformed equations of motion are 
integrated until steady state is attained. 

Most methods of generating systems of coordinates used in 
numerical solutions have been developed for elliptic problems. In 
these methods, the physical domain boundaries are known and the 
coordinate mesh is determined initially. Generally, the geometry 
of the mesh is not changed during the computation. Probably the 
most well known of these methods is the one developed by Thompson 
et al. (1) in which the transformed coordinates are obtained as a 
solution of Laplace’s equation in physical space. A number of 
other investigators (2, 3, 4) have developed schemes which can be 
used to generate appropriate coordinate systems using the same 
general idea. 
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Unfortunately, the solution of a separate elliptic equation 
is not conveniently included in the solution of a time-dependent 
set of equations. Hindman et al. (5) solved the two-dimensional time- 
dependent Euler equations with a truly adaptive grid scheme. The 
grid motion in time was generated by taking the time derivative of 
the governing differential equations of the coordinate mapping which 
was the same as that developed by Thompson. This provided the neces- 
sary grid speed equations which were then integrated to obtain the 
grid motion as a function of time. Hindman's work did not consider 
techniques which might be used to modify the location of the interior 
points depending upon the local solution. The interior point motion 
depended solely upon boundary motion. 

A technique for locating mesh points according to local flow 
information was presented by Dwyer et al. (6). This technique is 
similar to that used by Olson (7) and involves redistributing the 
mesh points at the end of any number of integration steps. This 
method does not permit a simple time integration of a differential 
equation similar to the equations of gas dynamics for the motion of 
the mesh points. It is the purpose of this paper to introduce a 
new technique which provides a simple way of moving the mesh points 
in physical space and reduces the error in the solution relative to 
that obtained using a fixed mesh. 

Pierson et al. (8) have also worked on the generation of grids 
which minimize error, but their technique involves the solution of 
a minimization problem. The extension of such a method to higher 
dimensions with the accompanying increase In the number of mesh points 



is not feasible due to the large amounts of computer time necessary 
to solve minimization problems. The method to be discussed in this 
paper is very simple in application and takes only a fraction of the 
time necessary to solve a minimization problem. 


THE METHOD 


To describe the basic idea employed in this paper, we consider 
transient problems in one space dimension. Let the physical space 
coordinates be x and t and let the computational space coordinates 
be £ and T where 
T = t 
£ = £(x,t) 


We require the calculation of the absolute value of the deriv- 
ative (ju^J) of some representative physical quantity (u) such as 
velocity, pressure, or temperature and the average value of the same 

derivative (lu r | ) for all mesh points. Given a certain number of 
1 £ 1 av 

grid points, truncation error can be minimized by allocating a number 

of points to the regions of large gradients and fewer points to the 

regions of small gradients. For an equispaced grid, a relocation of 

points in order to minimize error can be carried out. This can be 

achieved if points at which |u^| is larger than |u^|^ attract 

other points and points at which |u r | is smaller than |u~| repel 

s s av 

other points. In other words, every point induces a velocity at 
every other point, the magnitude and direction depending upon the 
local 'excess gradient'. It is logical to assume that the further a 


point A is from a point B, the smaller the effect of point A on B. 


This suggests that a 1/r law should be used. From the above con- 
siderations, it is possible to write 
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where i is the point at which the velocity is being determined, 

T N T is the total number of grid points, is the distance between 

points i and j in (£,t) space and f K T and ! n T are constants. 

The value of K can be determined if the maximum velocity that any 

point can achieve is specified. Convergence of the grid to a steady- 

state configuration is obtained by specifying a maximum value for 

K (K ). 
max 

Strong analogies can be found between the present formulation 
and treating the grid points as point electrical charges whose indi- 
vidual charges are proportional to the local T excess gradient 1 . 

The charges move so as to minimize the quantity 


: ° § [ ,u £ , j - KU] 


the minimum value of E being zero. 
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The collapsing of two computational space points into one physical 
space point is not possible because of two reasons: 

(a) The driving force g. 


8 “ lu d ' K'av (3) 

becomes negative when two points get very close and, hence, the points 
begin to repel each other. 

(b) The term ? x in Equation (2) gets very large as two points 
get very close. Hence, for a finite (£ ) T , tenc ^ s t0 zero ; i.e., 

the closer two points get to each other, the more difficult it becomes 
for them to move toward each other. However, Equation (2) does not 
prevent extreme stretching of the mesh in physical space, thus giving 
rise to errors in the calculation of the transformation metrics. The 
details of preventing extreme stretching for the problems solved in 
this paper are presented in the section on results. 

In the above discussion the driving force g is defined in terms 
of local and average first derivatives. A better formulation would be 
one in which g is defined in terms of quantities which are more 
representative of truncation error. One such quantity is the third 
derivative of u instead of the first derivative. The appropriate 
choice depends upon the order of the method being used and the problem 
itself. The flexibility in choosing the driving force and the quan- 
tity to be minimized is a particularly attractive feature of the 
current scheme. 

Two constants K and n appear in Equation (1) and a third one, 

K defines the maximum value that K can assume. The constants 
^ max 
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K and K together determine the grid speed. When K is less than 
max 

K , the grid speed is determined by K alone and when K is greater 
max 

than K , the grid speed is determined only by 1C . At present these 
max ° max 

constants are chosen empirically. In choosing these constants one 
should bear in mind that very large values of K result in grid 
oscillations which in turn result in longer convergence times, and very 
small values of K result in low grid speeds and hence, once again 
longer convergence times are observed. The constant K is calculated 
by knowing the maximum velocity that any point can achieve in the com- 


putational space 


i«i>j 


max 


The rules that govern the choice of 


[ (£.) 1 are the same as those that govern the choice of K 

max & max 

A variation of the constant f n* between 1 and 8 did not make any 
difference in the final grid in the one-dimensional case studied and a 
small difference in the two-dimensional case. The number of itera- 
tions for convergence increases slightly when larger values of n are 
used. However, larger values of n imply a smaller range of influence 
for any given point. Consider a value of n, 


n - 


log (2) 


When r = 2, 
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This implies that only points adjacent to a given point make a signif- 
icant contribution to the velocity of that point. Hence, Equation (1) 
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becomes 
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(4) 


The use of Equation (4) instead of Equation (1) greatly speeds up the 
grid generation process. 


EXTENSION TO MULTIDIMENSIONAL PROBLEMS 


The method can be extended to problems in two and three space 
dimensions without any difficulty. In particular, for a problem in 
two space dimensions, let the physical coordinates be given by (x,y,t) 
and the computational coordinates by (£,n,T) where 
T = t 

£ - £(x,y,t) 


0 = n(x,y,t) 


We now require the calculation |u^j and |u^| for every point and 


£ av 


for every row of points and u 
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for every column of points 


as in Figure 1. The grid speed equations are given by 
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r = \i (i-k) 2 + (j-A) 2 


where K^,!^ and n are constants, N the number of points in the 
E, direction and M the number of points in the rj direction. The 


values of and ^ can 


be determined by specifying (£. . ) 

I i j J ~ ■ 


max 


and | (r|^ j) T j nax respectively. Grid convergence can be achieved by 

specifying (K, ) and (K„) as in the one-dimensional case. 

r 1 max 2 max 

We also have the relationships 


(S. .) = (C x + C y ). . 

T VS X T Vt X,J 


(n. . ) = (n x + n y) . . 

1,J T XT y T 


which yield 


(6) 


(x ). . = 

ri,j 
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From Equation (7) it can be seen that the collapsing of mesh points 
and the overlapping of grid lines is again prevented as in the one- 
dimensional case. 

Points lying along a constant r| line can be made to move tan- 
gential to this line by specifying (q . .) to be zero for all these 

1,3 T 

points. A similar procedure can be adopted for constant £ lines. 
This facilitates the movement of points along surface boundaries, etc. 
However, this type of unnatural constraint on the velocity of points 
leads to a slightly distorted grid as shown in Figure 2. A more 
natural way of making points move tangential to boundaries is to 
specify periodic boundaries and use the pseudo points outside the 
region of interest also to calculate the grid speed. This procedure 
of calculating the grid speed results in the grid shown in Figure 3. 
The distortions present in Figure 2 are absent in Figure 3 and the 
grid is seen to be smooth and uniform. The grids shown in Figures 2 
and 3 were generated using a known solution to the two-dimensional 
transient, linear, viscous Burger's equation. 
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RESULTS 




The first problem solved using the present grid generation 
technique was the one-dimensional unsteady viscous Burger's equation 


u t + uu x = 


with the initial condition 


u(0,x) 


1 x = 0 

0 0 < x < 1 


and the boundary conditions 


u(t,Q) = 1 

u(t,l) = 3 

This problem has the steady state solution 

f Re 1 

u = u tanh — (1-x) 

l 2 J 


( 8 ) 


( 9 ) 


( 10 ) 
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where 


Re ■ 1/JJ 


ORIGINAL PA' R IS 
OF POOR QUALITY 

( 12 ) 


and u is the solution of the equation 

3-^ - exp {-uRe} ( 13 ) 

u+1 

The slope of the steady state solution at the right end increases and 
that at the left end tends to zero as Re increases. 

\ 
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McCormack’s method was used to integrate Equation (8) and three 
point central differences were used to calculate the metrics of the 
transformation. The stability limit for McCormack’s method for this 
problem was determined using the empirical formula given by Tannehill 
et al. (9). 

Results are presented for various values of Re in Figures 4-8. 

In all cases the steady state results using an adaptive grid and those 
obtained using an equispaced grid are compared with the exact solution. 
In Figure 4 results for Re = 1 are shown. The errors are very 
small ( < 0.04%) in both cases but the peak error without an adaptive 
grid is about 1.82 times the peak error with an adaptive grid. In 
Figure 5 results are presented for Re * 2. The ratio of the peak 
errors is now about 4.90 and a significant improvement in accuracy is 
seen. However, In Figure 5, the adaptive grid shows a slightly larger 
error in the region 0 < x < 0.2. This is due to the fact that the 
second point in the grj.d has moved to the right a substantial distance 
resulting in a higher error in this region. 

Figure 6 presents results for Re - 3. The second point in this 
case moves so far to the right that the truncation error In calculating 
the transformation metrics in this region swamps the entire solution 
resulting in a solution that is worse than the one obtained using an 
equispaced grid. In order to prevent extreme stretching of the grid 
it is necessary to include a measure of the truncation error Intro- 
duced in calculating the transformation metrics into the driving 


force g, 



g = [x^j - I x r | + 0 { |u r j - ju r | } 

fa t, 1 £ 1 av 1 5 1 £ 1 av 


( 14 ) 


where 0 is a constant. Since is greater than zero and is 

less than zero in this problem, Equation (14) can be written as 


8 ■ x ; - ( Vav • 8 {u | - (u e ) av > 


(15) 


Since the grid converges when g is a constant over the entire region, 
the transformation for the converged grid can be shown to be 


1 - fu - (l-f)(l-x) 0< f < 1 


(16) 


where f is a constant. Hence, an equivalent way of preventing ex- 
treme stretching is to define u as 


u = fu + Cl— f ) (1-x) 


(17) 


and the driving force g as 


£ av 


(18) 


The error curve obtained for Re = 3 and f = 0.7 is also shown in 
Figure 6, A substantial decrease in error is seen, the ratio of the 
peak errors being about 3.80. Figures 7 and 8 present results for 
Re “ 5 and Re - 10 respectively. In both cases a smoothed form of 
the solution as given by Equation (17) is used. The ratio of peak 
errors is about 2.23 for Re = 5 and 2.13 for Re * 10. Figure 9 
shows the transformation obtained for the case Re = 3, f = 0.7. 

The uniform nature of the transformation is apparent. 
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A better measure of the total truncation error at a point (e) is 


e a dx u 

XXX 

which can be approximated in this case as 
. 2 

e « dx u 

X- 

which yields 

e “ u JK 
A x 

Equation (21) suggests a driving force of the form 


( 19 ) 


( 20 ) 


( 21 ) 


g = 




tu 5 /5 xlav 


( 22 ) 


Results of using such a driving force for the case Re = 3 are pre- 
sented in Figure 10. The errors obtained are comparable to the ones 
obtained using an optimal f. However, the advantage in using this 
new form of the driving force lies in eliminating the empiricism 
required in determining the optimal f. Similar results were obtained 
for all Re <5.0. Excessive stretching was once again observed for 
higher values of Re, indicating the inaccuracy in estimating the 
error. The analysis and results presented in this and the preceding 
paragraph show that the method is limited only by the accuracy with 
which the total truncation error at a point can be estimated. 

The second problem solved was the two-dimensional unsteady, 
linearized, viscous Burger* s equation 


u + u + u 
t x y 


M (u 
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XX 
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( 23 ) 


421 



in a square domain with the initial conditions 

1 - exp(Re(x-l)) 

u(x,0,0) = 14 * — 

(1 - exp (-Re) ) 

1 - exp (Re (y— 1 ) ) 

u(0,y,0) = 1 + 

(1 - exp (-Re)) 


u * 1 otherwise 


where 

Re = 1/p 


and the boundary conditions 

1- exp(Re(x-l)) 

u(x,0,t) = 1 + 

(1 - exp (-Re)) 

1 - exp (Re (y-1)) 

u (0 ,y , t) = 1 + 

(1 - exp (-Re)) 


u(x,l,t) = 1 


u(l,y,t) = 1 




( 24 ) 


(25) 


(26) 



This problem has the steady state solution 

(1 - exp (Re (x-1 ) ) ) (1 - exp (Re (y-1) ) ) 


u = 1 + 


(27) 


(1 - exp (-Re)) 


McCormack's method was used to integrate Equation (23) and three 
point central differences were used to calculate the metrics of the 
transformation. To prevent excessive stretching of the grid a smoothed 
version of the solution (u) 

s 
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u = £u + (1-f) (4-x-y)/2 


0 < f < 1 


(28) 


Is used to calculate the driving force. 

Figure 11 shows the grid obtained for Re = 5 and f = 0.3. 

The error is calculated at the points shown in Figure 11 and a linear 
interpolation is used to calculate the error at the points correspon- 
ding to the equispaced grid. The results are presented in Figures 12- 
15, at each y station. The adaptive grid yields slightly higher 
errors in the low gradient region as in Figure 12 and gradually pro- 
gresses to much lower errors in the high gradient regions as in 
Figure 15. The increases in accuracy are not as high as in the one- 
dimensional case, the main reason being the inaccuracy in establishing 
the local truncation error. One complication that exists only in two- 
and three-dimensional problems is the appearance of cross derivative 
terms in any estimate of the local truncation error. The absence of 
cross derivative terms in the present formulation of the grid genera- 
tion scheme is felt particularly at the point x = 0.8, y — 0.2 in 
Figure 12. This point has a large value of and a small value of 

Uy resulting in mesh clustering only in the x direction. However, 

the terms u and u are by no means small and hence due to 

xyy xxy 

large Ay in this region give rise to large errors. Future work 
with two-dimensional problems will require that the influence of cross 
derivative terms be included in the generation of grids. 
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TIME REQUIREMENTS 


The number of integration steps required for convergence is always 
greater with an adaptive grid because of the lower values of maximum 
allowable time steps associated with mesh clustering. The ratio of the 
number of steps required with and without an adaptive grid goes all the 
way from 3.4 for Re = 10 to 1.4 for Re 2 1 in the one-dimen- 
sional case and takes on a value of 2.3 in the two-dimensional case. 
However, time estimates will be given only on a per integration step 
basis. In the one-dimensional case the generation of the grid and re- 
calculation of the transformation metrics takes less than 10% of the 
time taken for integration. In the two-dimensional case, the genera- 
tion of the grid takes 25% and recalculation of metrics takes 70% 
of the time taken for integration. One of the reasons for the 
excessive time taken for the calculation of metrics is the presence 

of second derivatives like E , £ , r) , and T) , all of which need 

xx yy xx yy 

to be determined numerically. The absence of these second derivatives 
greatly speeds up the calculation of metrics. Furthermore, if the 
problem requires the recalculation of metrics even without an adaptive 
grid, as in shock fitting programs, the time required to use an adap- 
tive grid becomes very attractive. It must also be remembered that 
the percent extra time In this case is high because the equation being 
solved is very simple. Since the time for grid generation remains 
about the same in far more complicated problems, the present extra 
time for grid generation will be much less for such problems. 
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In conclusion, the major contributions of this paper are: 

(a) Formulation of simple first order partial differential 
equations for the grid point velocity in transient problems. 

(b) Significant error reductions for solutions of Burger's 
equation in one and two dimensions. 

(c) The use of local flow information and boundary motion in 
determining the interior grid point motion. 
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Fig. 5. Comparison of errors 
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r v g . 7. Comparison of error 
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Fig. 9. Converged grid for the one 
dimensional Burger’s equa- 
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Fig. 11, Converged grid for the 
two-dimensional Burger 
equation . 








Fig. 12. Comparison of errors, 
Re = 5.0, Y = 0.2. 



Fig. 14. Comparison of errors 
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